Estimation of the bending rigidity and spontaneous curvature of fluid membranes in 

simulations 
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Several numerical methods for measuring the bending rigidity and the spontaneous curvature of 
fluid membranes are studied using two types of meshless membrane models. The bending rigidity 
is estimated from the thermal undulations of planar and tubular membranes and the axial force of 
tubular membranes. We found a large dependence of its estimate value from the thermal undulation 
analysis on the upper-cutofF frequency gcut of the least squares fit. The inverse power-spectrum fit 
with an extrapolation to gcut — >■ yields the smallest estimation error among the investigated 
methods. The spontaneous curvature is estimated from the axial force of tubular membranes and 
the average curvature of bent membrane strips. The results of these methods show good agreement 
with each other. 

PACS numbers: 87.16.D-, 87.17.Aa, 82.70.Uv 



I. INTRODUCTION 

When amphiphilic molecules are dissolved into aque- 
ous environments, these molecules self-assemble into sev- 
eral types of characteristic structures such as spherical 
or cylindrical micelles, bilayers, and inverted hexagonal 
structures [ij-Q- Among them, bilayer membrane is the 
basic structure of cells and organella. In living cells, 
biomembranes are not only static walls that separate 
components but also dynamical objects playing functions 
such as the vesicle transport of proteins via membrane 
fusion and fission. 

On a micrometer scale, the lipid membranes can be 
considered as a continuum surface, where the membrane 
thickness can be neglected. The curvature free energy of 
a curved membrane is given by Q 



-{Ci + c2-Cor + KC1C2 



dA, 



(1) 



where Ci and C2 are the principal curvatures at each po- 
sition of the membrane. The coefficients k and R are the 
bending rigidity and saddle-splay modulus, respectively. 
The spontaneous curvature Co vanishes when lipids sym- 
metrically distribute in both leafiets of the bilayer. The 
last term in Eq. ([1]) is constant for a fixed topology 
(Gauss-Bonnet theorem). The bending rigidity and spon- 
taneous curvature are basic quantities to understand the 
membrane properties. In this paper, we study the nu- 
merical measurement methods of the bending rigidity k 
and the spontaneous curvature Co in simulations. 

Several methods have been used to measure the bend- 
ing rigidity k in experiments and simulations. They are 
classified to two groups: (i) One utilizes the thermal fluc- 
tuations of the undulation modes of the membrane sur- 
face. In experiments, the surface fluctuations are mea- 
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sured by light microscopy with vesicles, cells, etc 0- 
[lOj . Theoretically, the fluctuation S2ectrum is derived by 
the perturbations from planar [l|, [Tl|, [T^ , spherical [ij- 
[l6j . and cylindrical [l7l - [l9| membranes. In simulations, 
the fluctuation spectrum of planar membranes is widely 
used to measure k [12|, [20l - [2a | . The fluctuations of tubu- 
lar membranes have not been simulated as yet, whereas 
those of quasi-spherical vesicles are calculated in Refs. 
[l^ll^ll^]. (ii) The other utilizes force measurements. 
A tubular (tether) membrane is formed from a liposome 
by a mechanical force (induced by optical tweezers, etc.). 
The bending rigidity can be measured using the force 
strength and the surface tension of the vesicle [28l - [3^ . 
The stability and the shapes of tubular membranes have 
been intensively studied in theories [TtHiQ . 33 37| . In 
simulations, Harmandaris et al. (38| measured k from 
the axial force and radius of tubular membranes. Re- 
cently, K was also measured from the surface tension of 
the buckled membranes in simulations (39| . 

In living cells, biomembranes have asymmetric lipid 
distribution in two leaflets [1^ . Such asymmetry of the 
membranes yields a non-zero spontaneous curvature Co. 
For a closed membrane i.e. vesicle, a low flip-flop rate 
between the leaflets can result in an effective spontaneous 
curvature (TBI . The vesicle morphology is varied with Co 
[4l| - |43j . Recent experiments show that the spontaneous 
curvature Co is also induced by grafting polymers, by 
absorption of protein onto the membrane surface, or by 
other means [44| - i48| . Since many proteins were found to 
control membrane curvatures, much attention has been 
paid to the effects of the spontaneous curvature Co. To 
simulate the effects of Co , it is important to establish nu- 
merical methods to measure Co- In previous studies, the 
spontaneous curvature is estimated from the comparison 
of membrane shapes with the results of the continuum 
theory [il], [i^. In this paper, we propose two direct 
methods for measuring the spontaneous curvature. 

Many types of membrane models have been devel- 
oped for simulations from the atomistic scale to a large 
micrometer scale (see review articles [sol - ls^ ). Among 
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them, particle-based meshless membrane models [261154 
|6]| are suitable for studying the large-scale membrane 
dynamics including topological changes such as mem- 
brane rupture, fusion, and fission. In these meshless 
models, a membrane particle does not represent a lipid 
molecule; it represents a membrane patch consisting of 
many molecules. The membrane particles self-assemble 
to form vesicles and planar membranes owing to their 
attractive interactions. 

We employ two types of meshless models in this pa- 
per: the mis membrane model [ssl . [56j and a new spin 
membrane model. In contrast to the mis model, the spin 
model allows a finite spontaneous curvature Cq similar 
to Yuan's meshless model [6l|. In our meshelss models, 
the bending rigidity k and the line tension F of the mem- 
brane edge can be varied separately for wide ranges of 
the fluid phase. 

In Sec. UTJ the membrane models and the simulation 
methods are described. The measurements of the bend- 
ing rigidity n from the undulation mode analysis of pla- 
nar membranes is explained in Sec. Iml The stretching 
force measurement and thermal undulations of tubular 
membranes are described in Sec. IIVI The measurement 
of the spontaneous curvature Co from the force measure- 
ment of tubular membranes and the average curvature of 
bent membrane strips is explained in Sec. [V] A summary 
is provided in Sec. IVIl 



where is the thermal energy. The potential U con- 
sists of a repulsive soft-core potential J7repj an attractive 
potential C/att, with a coefficient e, and a curvature po- 
tential Ua. In a quasi-two-dimensional membrane sur- 
face, the particles interact via the potentials Uicp and 
t/att- The repulsive excluded interaction potential of a 
diameter a is given by 

C/rcp = ^exp(-20(r,j/a - 1) + B) f,^t{n, / a) , (3) 

i<j 

where B = 0.126, and rij is the distance between parti- 
cles i and j. The diameter a is employed as the length 
unit to display the simulation results throughout this pa- 
per. A C°°-cutoff function 

/cut(.) = ( .^"P^'^^' + W^)} < ^^-) (4) 
L t) [S ^ Scutj 

is employed. For Eq. ([3|), the values n = 12, a = 1, and 
Scut = 1.2 are used. 

A solvent-free membrane model requires an attractive 
interaction mimicking the "hydrophobic" repulsion be- 
tween hydrocarbon chains of lipid or surfactant molecules 
and aqueous solvent. We employ a potential 

C/att = J2 0-25 ln[l + cxp{-4(p, - p*)}] - C, (5) 



II. SIMULATION MODEL AND METHOD 

We employ two types of meshless membrane models. 
They use different curvature potentials, (i) The mesh- 
less mis membrane model (ssl . [56j : Membrane parti- 
cles possess only translational degrees of freedom, and 
form quasi-two-dimensional structures stabilized by a 
mutibody potential based on moving least-squares (mis) 
method jH,!!^. (ii) The meshless spin membrane model: 
Particles also possess orientational degrees of freedom 
and interact with potentials similar to those described in 
Ref. [1^. In both the models, the bending rigidity k and 
the line tension F of the membrane edge are controlled 
by changing the parameters of particle interactions (see 
Appendix). Further, the spontaneous curvature Co can 
be varied in the spin model. 



A. Meshless mis membrane model 

Since the model is explained in detail in Refs. [H^, 
[56j . the model is outlined only briefly in this section. 
A membrane consists of N particles, which possess no 
internal degrees of freedom. The particles interact with 
each other via a potential 



U 



e(C/rep + C/att) + C/a, 



(2) 



which is a function of the local density of the particles 
Pi, which are defined by 



(6) 



For this attractive interaction, the values Shaif =1.8 (at 
which /cut(shait) = 0.5), Scut = Shaif + 0.3, and n = 
12 are used. The factor a in /cut(s) is given by a = 

In(2){(seut/Shalf)" - 1} ^ 3.715. 

The constant C = 0.25 ln{l + exp(4p*)} ~ p* is set to 
achieve C/att = at pi = 0. For pi < p*, the potential is 
approximately C/att — —Pi, and therefore, it acts as a pair 
potential with C/att ^ - J2i<j '^fcutinj /a). For Pi > p* , 
this function saturates to the constant —C. Thus, it is 
a pairwise potential with a smooth cutoff at the density 
Pi = p* . Wc set p* = 6 in this paper to simulate a fluid 
membrane. 

In addition to the potentials C/rop and C/att, we add a 
curvature potential 



Ua = ka'^ap\{ri), 



(7) 



where the shape parameter "aplanarity" is defined by 
QD^ 9A1A2A3 



CKpl 



TwMw (Ai + A2 + A3)(AiA2 + A2A3 + A3A1) ■ 



(8) 

The aplanarity represents the degree of deviation 
from the planar shape, and Ai < A2 < A3 are the three 
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eigenvalues of the weighted gyration tensor for the zth 
particle given by a^p = J2ji^j - "G)(/3i - ^G)wcv{nj), 
where a, /? € {x, y, z} and the local center of mass 
= J2j ^jWcvinj)/ Wcvin,])- The aplanarity is cal- 
culated from three rotational invariants of the gyration 
tensor: the determinant D„, trace 
its three minors, = 



A,yy 



W 1 



and the sum of 



^yz~C'zx- This aplanarity ttpi takes a value in the interval 
[0, 1] and is proportional to Ai for Ai <C A2, A3. 

A Gaussian function with a C°° cutoff [55[ is employed 
for calculation of the weight of the gyration tensor, 



Wcv(r-,j) 



which is smoothly cut off at 




< rcc) 



(9) 



rameters n = 12, 



0.5rcc, and Vr. 



We use the pa- 
: 3cr here. When 



the ith particle has two or less particles within the cut- 
off distance rij < rcc, the particles could be on a plane, 
thereby making api = 0. 



B. Meshless spin membrane model 

In the meshless spin membrane model, a tilt potential 
Unit and a bending potential /7bend are employed instead 
of the aplanarity potential Ua- U/k-QT = e(J7rop + C^att) + 
^^tiit + C^iend- In this model, the spontaneous curvature 
of the membrane can be varied. Each membrane parti- 
cle has one orientation unit vector itj {\ui\ = 1). The 
potentials Utm and J/bond are given by 



end 



^(Mi - Uj - Chdfi/fwcv{rij), (10) 



where = rij /rij. These potentials are the discretized 
versions of the tilt and the bending potentials of the tilt 
model (62l . l63j , respectively. The spontaneous curvature 
of the membrane can be controlled by the potential pa- 
rameter Cbd as discussed later in Sec. |Vl Recently, we 
employed similar potentials for a molecular lipid model 
to form bilayer membranes psj . In the molecular model, 
the positions = r, -I- Uicr were used for calculation of 
the weight Wcviftj) t° stabilize the bilayer structure. In 
contrast, here, the distance between the centers of 
mass of the particles is used to calculate Wcv{rij)- We 
use the same parameters in the functions C/rep, C^att, and 
Wcvijij) for both the models. Unless otherwise specified, 
we use e = 4. 



C. Simulation methods: Brownian dynamics 

We simulate the membranes in the NVT ensemble (the 
particle number, volume, and temperature are kept con- 



stant). The dynamics of the membrane is calculated us- 
ing Brownian dynamics (underdamped Langevin equa- 
tion). The motion of the membrane particles is given 
by 



(Pvi ^ dri r , N c^C/ 



dt2 



dt 



I- 



dri 

( , dU\^ ^ 



(11) 



12) 



where m and / arc the mass and the moment of in- 
ertia of the molecule, respectively. The angular veloc- 
ity u)i = dui/dt is rotated by the perpendicular force 
f ^ = f — (f • Ui)ui. The length u| = 1 is kept constant 
by a Lagrange multiplier A. The friction coefficients Co 
and and the Gaussian white noises gp(t) and gl{t) 
obey the fiuctuation-dissipation theorem; 



(13) 



(5t,(*)5,^o.(i')> = 2kBTCpAjSc.,o.Jp,pAt - t'). 

Here, ai,a2 G {x^y^z} and /3i,/32 {G,r}. In the fol- 
lowing sections, we use the time unit t = C^Qa'^ /k^T and 
the energy unit k^T . We use m = C,gt and / = Ci.r. For 
the mis or the spin model, Eq. pT|) or Eqs. (jlip and 
(fT2|) are integrated by the leapfrog algorithm with a time 
step of Ai ~ 0.005t, respectively. The simulations are 
performed with periodic boundary conditions in a box of 
dimensions Lx x Ly x Lz- 



III. THERMAL UNDULATIONS OF PLANAR 
MEMBRANES 



The undulation spectrum analysis of a planar mem- 
brane is the most widely used method to estimate the 
bending rigidity k in simulations. In this section, we com- 
pare fitting methods and establish a large dependence 
of the estimate value on the cutoff frequency (/cut- The 
spectrum of the undulation modes (|/i(g)p) of the planar 
membranes in a Fourier space is given by [l|, [l^ [20| 



kBT 



79 



(14) 



We calculate |/i((?)P for the planar membranes with 
Lx = Ly from the raw data (the particle position r^) as 
well as from the positions averaged on a \/N/2 x y/ N/2 
square mesh with {x,-nh,ymh) = (dmhn-s, dmh^^y)- The 
height Zmh of a mesh point is obtained from the weighted 
average of the molecular position in the four neighbor 
cells with z,nh = J2iZiWnihixi,yi)/{J2iWinhixi,yi)) and 

Wmh{Xi,yt) = (1 - \Xt - a;mh|/rfmh)(l - \yi - 2/mh|/rfmh). 

We refer to the former and the latter spectra as the "raw 
spectrum" and the "mesh spectrum" respectively. Fig- 
ure [TJa) clearly shows the dependence of the tension- 
less membrane (surface tension 7 = 0). For increasing 
7, (|ft-((7)p) decreases at a low frequency q (compare the 
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FIG. 1: (Color online) Spectra of undulation modes 
of a planar membranes for the mis membrane model with 
ka = 10 and e = 4. Results for calculated from the 

particle positions (raw; +) and from the averaged positions 
on a square mesh (mesh: x ) are shown, (a) Spectra for 7 = 
{A^ylNa'^ = 1.416) and 7 = 3 {A^.yjNa'^ = 1.55) at the 
number of particles A'' = 6400. The thick black line represents 
a slope of The inset shows (|ft(g)|^) of raw data for a 

large q at iV = 400, 1600, and 6400. The dependence of 
l/(|/i(g)P) on for iV = 6400 and TV = 1600 is shown in (b) 
and (c), respectively. 



data at 7 = and 7 = 3 in Fig. HJa)). The mesh spec- 
trum approaches zero at a high q. On the contrary, the 
raw spectrum saturates at a finite value, which increases 
with increasing N (see the inset of Fig. [TJa)). These 
high q modes are caused by particle protrusion due to 
the short-range potential interactions between particles. 
Averaging over the mesh removes most of the effects of 
these particle protrusions at a high q. The effects of pro- 
trusions on bending rigidity estimation will be discussed 
in the last part of this section, with mis and spin models. 

For the estimation of k, two types of fits are widely 
used, (i) a log-log fit for the tensionless membranes. 



ln((IM9)P)) 



In 



/ K 



VfcBT 



41n((z), 



and (ii) an inverse power-spectrum fit, 
1 



79 



nq 



(15) 



(16) 



In the latter fit, 7 can also be estimated. We calculated 
both of them using the linear least squares fit for vari- 
ous cutoff values given by gcut- For the fits with one fit 




qcuA 



FIG. 2: (Color online) Estimation of the bending rigidity /t of 
the tensionless membrane from the log-log fit for (a) the mis 
membrane model with ka = 10 and (b) the phenomenological 
function Eq. p9l) . (a) The symbols (blue or dark gray: □, V) 
and (red or light gray: o, A) represent the data for A'^ — 
1600 and 6400, respectively. The symbols at relatively lower 
positions (V, A) and the symbols at relatively higher positions 
(□, o) represent the data estimated from the fits for the raw 
spectrum and the mesh spectrum, respectively. 



parameter (k), the surface tension calculated from the 
pressure tensor is used as the value of 7. The surface 
tension is given by 

7 = {Pzz - {Pxx + Pyy)/2)L,, (17) 

with the diagonal components of the pressure tensor 



dU. 



Pec = (A^fcBT-^a,— 



(18) 



where a £ {x,y^z} [64l |65| . In calculating Paa, the peri- 
odic image ai + nLa nearest to the other interacting par- 
ticles is employed, when the potential interaction crosses 
the periodic boundary. We compare these two surface 
tensions, calculated from the pressure and the thermal 
undulations, later in this section. 

Figure [2ja) shows the bending rigidity k estimated 
from the log-log fit of Eq. for the data for which 

Q < ^cut- As the cutoff frequency gcut increases, the es- 
timate value of K gradually decreases. The fits to the 
mesh spectrum are less sensitive to gcut than those to 
the raw spectrum. This gcut dependence is caused by 
the neglected fluctuation modes. As shown in Fig. [TJa), 
(|ft.(g)P) of the raw spectrum deviates from g^'* at a high 
g because of particle protrusions. Goetz et al. reported 
that the protrusion modes of lipid molecules have a q~^ 
dependence [l3] . To clarify the influence of these protru- 
sion modes on the spectrum, we test the least squares fit 
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to a phenomcnological function 

ilhil)?) -—{r^^Y--, + (19) 

The first term is the thermal-undulation mode with the 
bending rigidity ko and surface tension 70 = 0, and the 
last term (oc q^^) is a protrusion mode. Figure [Hb) 
shows K estimated by the log-log fit of Eq. (|16p to the 
phenomcnological function of Eq. (fT9)) with 70 = for 
O.OItt < g < gcut- At a high 6prt or high gcut, the fit 
gives a lower k than the actual bending rigidity kq- It 
qualitatively reproduces the (/cut dependence obtained in 
the simulations. A similar decrease is observed when 
^prt/?^^ is used instead of bprt/q~^- Thus, this depen- 
dence is not sensitive to the function shape. As a result, 
the fit with a lower gcut should yield a more accurate k, 
close to kq. However, since the number of the data points 
for the fits is fewer for lower gcut and the spectrum at low 
q has a larger statistical error, the error bar is larger in 
low gcut region. Thus, the medium cutoff (qcut — O.Itt 
in Fig. [5]) is a reasonable choice for the estimation of k 
from the log- log fit. 

Next, we explain the inverse power-spectrum fit of Eq. 
([T6| . This fit also shows a large dependence on gcut- 
Four types of fits (with one or two fit parameters for the 
raw and mesh spectra) for the tensionless membranes 
are shown in Fig. [31 The spectra are fitted with one 
(k) or two (k, 7) fit parameters. The one-parameter fit 
provides a lower slope of the Qcut^f^ curve than the two- 
parameter fit. The mesh spectrum gives lower slopes for 
K and 7 than the raw spectrum does (see Figs, ^a) and 
(d)). Thus, substantial difi'erences are seen between the 
K values estimated by different fits for a finite Qcut- How- 
ever, all the fits converge at gcut — > 0. A similar gcut 
dependence is observed for the fit to the phenomcnolog- 
ical function in Eq. (fT9|) with 70 = 0. At gcut 0, k 
values converge to the correct value kq- Therefore, the 
bending rigidity n can be estimated by an extrapolation 
using the linear least squares fit to a straight line (see 
solid lines in Figs. EJa) and (b)). 

We select the extrapolated value at Qcut = from a 
one-parameter fit for the mesh spectrum as the bending 
rigidity k, considering that it has the lowest slope of the 
qcut-K curve. We estimate the numerical error in k from 
two contributions; i.e., for Ak = Ak{q -|- A^ti: the maxi- 
mum difference Ak{q between four extrapolated values is 
considered as the numerical error size from the choice of 
functions, and A^fi is the error of the least squares fit. 

As the membrane area increases, the surface tension in- 
creases. We investigated the area dependence of k using 
the extrapolation method for the parameters {k^ = 10, 
e = i) investigated in our previous paper [s^ . The intrin- 
sic area A of the membrane is larger than the projected 
area A^y in the xy plane because of the membrane un- 
dulations. We calculate A from the mesh points for the 
mesh spectrum. Figure H] shows the area dependence of n 
and the difference A7 = 7fl — 7pr in the surface tension es- 
timated from the two methods. The surface tensions 7pi. 




FIG. 3: (Color online) Estimation of the bending rigidity k 
from the fit to l/{\h{q)f) with Eq. ^ for (a) the mis model 
with ka ~ 10, (b) the spin model with fcbend = ktnt = 15 
and Cbd = at e = 4, and (c) the phenomenological function 
in Eq. (|19|l. The symbols (blue or dark gray: □, *, +, V and 
red or light gray: 0,0, x , A) represent the data for A*' = 1600 
and 6400, respectively. From the top to the bottom, the one- 
parameter fit for a mesh spectrum (njo); two-parameter fit 
for a mesh spectrum (*,o); one-parameter fit for a raw spec- 
trum (-I-, x); and two-parameter fit for a raw spectrum (V, A) 
are shown. The solid lines in (a) and (b) are obtained by 
a least squares fit for the data in (a) ((/cut/Tr)^ < 0.08 and 
(b) (gcut/7r)^ < 0.04. (c) The solid and dashed lines repre- 
sent the data for 70/Ko = and 0.05 at bprt = 1, respec- 
tively. The upper or lower lines show the results obtained 
using one- («) or two- (k, 7) parameter fits, respectively, (d) 
The surface tension estimated from the two-parameter fit for 
the mis model. The solid lines represent the curves fitted by 
7 = 7fl + a.^iqcMt + &7(7cut for the data at (gcut/vr)^ < 0.08. 
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FIG. 4: (Color online) Intrinsic area A dependence of (a) the 
bending rigidity k and (b) the difference A7 = 7fl — 7pr of the 
surface tensions estimated from the undulations (Eq. (|16[)) 
and the pressure tensor (Eq. (flTt '). The solid line with (□, V) 
and dashed line with (o, A) represent the data for A'^ = 1600 
and 6400, respectively, (b) The symbols (V,A) and (0,0) 
represent the data estimated from the raw spectrum and the 
mesh spectrum, respectively. 



are estimated from the pressure tensor (Eq. ^(T7\ ): jp^ = 
and A/Na^ = 1.443 at A^y/Na^ = 1.416, whereas 
7pr = 2.96 and A/Na^ = 1.565 at A^y/Na^ = 1.55. 
The surface tension estimated by the extrapolation with 
7 = 7fl + a7'7cut + ^79cut has a very good agreement with 
7pr, as shown in Fig. IH^b). Previously, Imparato reported 
that the surface tension estimated from the thermal un- 
dulations is smaller than that estimated from the pres- 
sure tensor in molecular simulations [66| . The reported 
difference may be due to the similar effects of a finite 

9cut- 

Although the estimated k decreases with an increase 
in A, it is accompanied by larger error bars. Thus, we 
do not find a clear dependence of k on 7. For a finite 
7o, the fit for the phenomenological function of Eq. 
shows a deviation from the correct value kq at gcut 
(see Fig. [3{c)). The one-parameter fit shows an abrupt 
decrease at (/cut — 0, whereas rapid decreases arc not 
observed in the simulations. These results suggest that 
systematic errors may be involved in the estimation of k 
at a finite 7. Further investigations are needed to clarify 
the K dependence on 7. 

The extrapolation method works well for the spin 
membrane model as well as the mis membrane model. 
The spin model has a larger (jcut dependence. The gcut^^ 
curve deviates from the straight line at (gcut/Ti")^ ^ 0.05 
and (gcut/7r)^ ^0.1 for the spin and mis models, re- 
spectively (compare Figs, ^a) and (b)). This suggests 




FIG. 5: (Color online) Snapshot of a tubular membrane in the 
simulation of the mis membrane model at = 2400, kc, = 10, 
e = 4, and Lz — 80a. 



that the spin model has greater particle protrusion than 
the mis model. Because the protrusion is induced by the 
short range interactions between particles or molecules, it 
is sensitive to the potentials of simulation models. To use 
this method, one should ensure that the spectrum at a 
sufficiently low q is included for the linear extrapolation. 
The bilayer membranes of the spin molecular model (25j 
show a similar dependence as that of the meshless spin 
model (data not shown). Simulation models accompa- 
nied with larger protrusions require larger system sizes 
to estimate k from the thermal undulations. 

For K extrapolated at gcut 0, no significant depen- 
dence on the system size N is detected. All of the ex- 
trapolations for = 1600 and 6400 converge (see Figs. 
[3]and|3]). At a finite q, the raw spectrum is dependent 
on N because the protrusion amplitude increases with 
increasing N (compare Figs. [T{b) and (c)). The slope 
of the qcut~K curve is higher at a larger N for the raw 
spectrum (see Fig. [3]). 

We also estimated k using a nonlinear least squares 
fit for {\h(q)\'^) without and with the protrusion term 
^ q-^ [Eq. (dH) and Eq. (HH)]. Although these fits 
are not sensitive to (jcut, they have larger errors than the 
above fitting methods. Therefore, we conclude that the 
inverse power spectrum with an extrapolation to q^ut ^ 
is the best fitting method for the undulation of planar 
membranes. In the Appendix, we list the bending rigidity 
K estimated by the inverse power-spectrum fit for the mis 
and spin membrane models with various parameters. 



IV. TUBULAR MEMBRANES WITH NO 
SPONTANEOUS CURVATURE 



In this section, we present the estimation of the bend- 
ing rigidity k from tubular membranes. For a tubular 
membrane with a radius R and a length Lz, the curva- 
ture free energy Eq. ([1]) is written as 



J" = 2ttRLz 



Cr 



(20) 
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Under the fixed area condition A = 2-kRLz 
axial force = dF/dLz\A is given by 



const, tlie 



(21) 



lateral tension is anisotropic: = /z/'^ttR = k{1/R — 
Co) in the axial direction, 79 = in the azimuth direc- 
tion, and 7av = 7^/2 in average. Although we assume 
the constant area here, the area compressibility does not 
change the force, Eq. ([21]) . When the area compressibil- 
ity is taken into account, the free energy Eq. ([20l) has 
an additional term C/ar(A) {Uar{A) = Ka{A - Ao)V2^o 
ioi A — Aq ^ 1, where Aq is the area of the tension- 
less membrane). Then, the same force is derived from 
fz = dF/dLz and dF/dR = 0. At Co = 0, the force is 
inversely proportional to i? (/^ = 271 k / R). Using this re- 
lation, K was previously estimated in experiments (28i - [3^ 
and molecular simulations [ssj . 

In this section, we investigate stretched cylindrical 
membranes with Co = using the mis membrane model 
(see Fig. [5]). All the tubes are connected periodically 
in the axial direction with the periodic length L^- The 
initial conformations at each Lz are made by slow stretch- 
ing or shrinkage of the length Lz with a speed less than 
dLz/dt = 0.002cr/r. We checked that no hystereses are 
seen in results between stretching and shrinkage. The 
bending rigidity k is measured at fixed Lz- After dis- 
carding the data for the first calculation period 1600t, 
the data are averaged for a period 48000t (72000t) for 
N = 2400 (1200). Eight simulations starting with inde- 
pendent initial conformations are performed. 

Figure|n]shows the estimate values of the bending rigid- 
ity Kest = fzR/'^T^ for various Lz- The inverse radius of 
the tube a/R is employed for the horizontal axis. With 
increasing cr/i?, the cylinder tube becomes narrower and 
longer. The radius R of the cylinder is simply estimated 
by averaging the distance of each particle from the cylin- 
drical axis: R = (E,;{(^,: - ^g)' + {y^ - yc^V^VN), 
where {xcUq) is the center of mass of the membrane 
projected on the xy plane. The estimate values of Kest 
for long tubes {(j/R > 0.06) have very good agreements 
with those obtained for planar membranes in Sec. IIIII 

For shorter tubes {a/R < 0.06), large bumps (or 
peaks) appear at = 1200 (see Fig. IH^a)). Because 
they are suppressed at A^ = 2400 (tubes that have twice 
the length for the same radius), these bumps are likely 
caused by the finite size effects in the z direction. For 
long tubes, Kest decreases slightly with increasing a/R. 
This may show the dependence of k on the area expansion 
or be on account of the higher-order terms of the bend- 
ing elasticity discussed in Ref. [s^. As explained later 
in this section, the intrinsic area A is larger for longer 
tubes. The Kest decrease resembles that seen in the esti- 
mation for the planar membranes of Fig. U) The decrease 
rate strongly depends on ka; —d{Kcst/k^T)/d{a/R) = 3, 
11, and 31 for = 5, 10, and 20, respectively. This 
difference can be partially explained by the effects of the 
thermal fluctuations. Recently, Barbetta et al. [l^ de- 
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FIG. 6: (Color online) (a) Bending rigidity /tost = fzR/'iir es- 
timated from the stretching force measurement using the mis 
membrane model. The symbols (blue or dark gray: -|-,V,o 
and red or light gray; o. A, x ) represent the data at A'^ = 1200 
and 2400, respectively, for ka = 20 (-l-,o); 10 (V, A); and 5 
(o, x). (b), (c), and (d) show magnified plots of the data 
shown in (a) for ka = 20, 10, and 5, respectively. 



rived the axial force under the thermal fluctuations using 
perturbation theory. In their theory, the force is given by 



2ti 



27r2 



(22) 



where A is the cutoff wave vector. For A = 1/2ct, Kest = 
!zR/2tx increases by k^T from a/R = 0.1 to 0.25. Thus, 
the decrease rate for fc„ = 5 in Fig. [Hlis reduced by the 
thermal-fluctuation effects. 
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7: (Color online) (a) Undulation spectra (|it„ 
the tubular membranes of the mis model at A'^ = 2400 and 
= 10, for LJa = 112 (A), 64 (v), 40 (x), and 20 (o). 
The data for < m, n < 10 are shown. The horizontal axis is 
the fourth power of the effective frequency Q defined in Eq. 
(I25j. (b)(c) Spectra for (b) m = and < n < 10 and (c) 
< m < 10 and n = are extracted from the data in (a) to 
clearly show the m and q dependence. 



The bending rigidity k can also be estimated from 
the surface undulations of tubular membranes. Recently, 
Fournier et al. studied the thermal undulations on a 
cylindrical membrane theoretically [l8| . They predicted 
the nontrivial effects of the critical Goldstone modes for 
narrow and long tubes. We numerically analyze the sur- 
face fluctuations and membrane area in comparison with 
their theoretical framework. The membrane position 
is expressed in the cylindrical coordinates r{d,C,)/R = 
([l + w(e',C)]cos6', [l + u{9,C)] sm9,C) for < 6* < 27r and 
< C < Lz/R- We calculate u{9,C) from the raw data 
of the particle positions r^. The cylindrical axis {xq, Vg) 
and the radius R are estimated in the same manner as in 
the above force measurement. The Fourier modes of the 
cylindrical surface fluctuations are given by 



<0,O = 



R 



2ttL, 



J{me+qC) 



(23) 
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FIG. 8: (Color online) Estimation of the bending rigidity k 
for the mis membrane model from the fit to l/(|Mm,q|'^) with 
Eq. for L^/a = 80 (□), 64 (v), 40 (x), and 32 (o), 

for various cutoff frequencies Qcut- The data are fitted for 
< m, n < 10 and Q"^ < Qi^^ 



where q 2TTnR/Lz, \m\ < 2ttR/1, and \n\ < L^/l. The 
cutoff length I is the mean distance between neighboring 
membrane particles in meshless membrane models or the 
membrane thickness in molecular models. 

At thermal equilibrium, surface undulation of the 
cylindrical membrane can be estimated by the pertur- 
bation theory. The spectrum is given by [l3, [l3 



{\Un 



2m^ 



(24) 
(25) 



where Q denote the normalized amplitude of frequen- 
cies in the two-dimensional cylindrical space. Since the 
expression Eq. (|24p . for cylindrical membranes, is the 
counterpart of Eq. (|14p for planar membranes, k can be 
estimated using a similar fitting method. 

Figure [7] shows the undulation spectra for L^/a = 20, 
40, 64, and 112 {R/a = 27.1,13.6,8.58, and 5.01) at 
N = 1200 and ka = 10. The modes for small frequencies 
at < m,n < 10 are shown here ((|wi_oP) are omit- 
ted because of their divergence). While the amplitudes 
(|um,„P) at a low q exhibit deviations for the narrow 
tubes [Lz/a = 112), they clearly show dependence 
at a low Q when the ratio between the circumferential 
length and the cylinder length ^-KRjLz is around unity. 
In these regions, k can be estimated by least squares fits, 
as explained in Sec. IIIIl 

Figure [8] shows the bending rigidity k obtained by a 
linear least squares fit with various upper cutoffs Q; 
the inverse power spectrum (|um,9p)"^ using Eq 
in a manner similar to that used in Sec. III. When the 
horizontal axis is normalized by {2'kR/ , the data for 
all Lz overlap. This dependence is very similar to that of 
the planar membranes shown in Fig. Elja). k approaches 
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FIG. 9: (Color online) Dependence of estimate values of the 
bending rigidity Kfiuc from the surface fluctuation spectrum 
on the frequency Q*, for L^/a = 112 (A), 64 (v), 40 (x), 
and 20 (o), at kc = 10 and = 2400. 
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FIG. 11: (Color online) Tube radius dependence of the intrin- 
sic membrane area A and the projected area = 2nRLz at 
ka ~ 10. The area is normalized by that of the tensionless 
planar membrane {Ao/Na^ = 1.443). Blue (dark gray) and 
red (light gray) points show the results for = 1200 and 
2400, respectively. 



^ I I I I I I I I I I I I I I I 

0.8 — — 

_g 

it ~ 

0.6 — — 

_l I I I I I I I I I I I I I I I I I I 

0.05 0.1 0.15 0.2 0.25 

a/R 

FIG. 10: (Color online) Ratio Kfluc/^^cst between the estimate 
values of the bending rigidity from the force measurement Eq. 
(|21[l and from the surface fluctuation spectrum shown in Fig. 
|9] at fco = 10 and iV = 2400. The lowest ten modes are 
used to estimate khuc- 

a value of around 20 as Qcut 0. Thus, the spectra 
of {\um,q\'^) well reflect the bending rigidity of the mem- 
branes. 

Figure O shows the dependence of estimate values of 
the bending rigidity kauc = / {\um.:q\'^)Q'^ on the fre- 
quency (Q^) using Eq. In this figure, we plot the 
data for < m,n < 10. In intermediate length scales, 
where the radius of the tube is 0.045 < a/R < 0.1, the 
spectrum of kauc collapses into a smooth shape, whereas 
systematic deviations for a specific to or g appear at 
Lz/<7 = 112 and 20. Low Q values of kauc, which rep- 
resent longer-wavelength properties, well converge to a 
value around Kost- To compare kauc with Kestj their ra- 
tio Kfluc/'^ost is shown in Fig. 1101 In the region plot- 
ted in Fig. [To] (15 < Lz/cr < 25), the ratio is constant 
~ 0.9, whereas both of them decrease with increasing 
a/R. These tendencies are also obtained for other ka or 
N. The results could be reflecting the dependence of k 
on the intrinsic area A for the meshless model. 

In the meshless membranes, the tubular membrane 
area is slightly expanded owing to the axial tension. Here, 



we estimate the membrane area in the following manner: 
A Delaunay tessellation is performed for the {9,C) coor- 
dinates to construct a triangulated surface on the mem- 
brane. Then, the intrinsic membrane area A is calculated 
as the sum over the area of the triangles for the 3D par- 
ticle positions. Figure [TT] shows the intrinsic area A for 
N = 1200 and 2400 in comparison with the projected 
area = 2tt RLz- As the membrane area is expanded 
for larger axial tension (7^ = k/2R^), A and approach 
each other. The area expansion of A is more than twice 
that of the planar membranes for the same average sur- 
face tension 7av. The anisotropy of the surface tension 
results in a low effective area compression modulus Ka- 
Fournier et al. (isj derived the dependence of the nor- 
malized excess area 

(26) 

on the axial tension c^, from the undulation spectrum 
Eqs. (HH) and ([25]), where A A = A - A^. They pre- 
dicted that a higher axial tension generates an increase in 
the normalized excess area owing to the enhanced Gold- 
stone mode fluctuations, contrary to the case for planar 
membranes. Figure [12] shows the Lz/R dependence of 
the excess area obtained from the perturbation the- 
ory [Eq. (|26|) ] and the simulations. When Ocx is plot- 
ted for Lz/R, the size dependence of Oox from Eq. 
is seen only for the middle region of Lz/R ^ 10, and 
then, all the three curves converge at Lz/R —5- and 
Lz/R ^ 00. The enhanced fluctuations in the azimuth 
or axial direction generates a large Oox at Lz/R ^ or 
Lz/R — >■ 00, respectively. For Lz/R < 10, our simulation 
shows good agreement with their prediction. However, 
for Lz/R ^ 10, aex decreases in the simulation but in- 
creases in their theory. Thus, the thermal undulations 
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FIG. 12: (Color online) Excess dependence on L^/R 

calculated from the simulations and the perturbation theory 
Eq. (|26|) . For the simulations, flex is shown in the range 
10 < L^/a < 72 and 18 < L^/a < 144 for TV = 1200 (o) and 
2400 (x) at ka — 10, respectively. The analytical data from 
Eq. (123 are also shown for L^/l = 10^ (+), 10^ (A), and 10^ 
(□). 

of longer tubes are suppressed in the simulations. This 
discrepancy may be caused by the suppression of the pro- 
trusion modes or the effects of the higher-order terms of 
the perturbations. Further study is necessary to clarify 
the origin of this difference. 

When the solvent is explicitly taken into account or 
when bilayer membranes have a low flip-flop frequency, 
the bending rigidity is difficult to simulate using tubular 
membranes. The tubular membranes would, in such a 
case, exhibit very slow relaxation to the thermal equi- 
librium state, considering that a radius variation of the 
tubular membrane accompanies changes in the tube vol- 
ume and in the area difference between the two leaflets. 
Therefore, the Laplace pressure needs to be taken into 
account or an additional numerical technique is required 
to exchange the solvent particles or lipids between the 
upper and the lower sides of the bilayers. 



V. MEASUREMENT OF SPONTANEOUS 
CURVATURE 

Next, we investigate the estimation method of the 
spontaneous curvature Cq using the spin membrane 
model. We estimate Cq from the axial force of tubular 
membranes and the shape of membrane strips. When the 
membranes have a nonzero spontaneous curvature, the 
membrane tube has the lowest free energy at 1/R = Co, 
where becomes zero (see Eqs. (PH)) and (HJ)). Fig- 
ure [13] shows that 1/R-fz lines move down for increasing 
Cbd- The finite-size effects discussed in Sec. IIVI for a 
small 1/R are also seen for a finite Cq (See Fig. [BJa)). 
The spontaneous curvature Cq and the bending rigid- 




a/R 

FIG. 13: (Color online) Force /z dependence on the radius 
R of the membrane tube for the spin membrane model at 
e = 5 and Cbd = 0, 0.1, and 0.2. (a) fcbcnd = fctiit = 20. 
(b) fcbcnd ~ fctiit = 5. The symbols (o, V,n) and (x,A,o) 
represent A'' = 1200 and = 2400, respectively. The solid 
lines are obtained by linear least squares fits. 



ity K were estimated by the linear least squares fit to 
Eq. dH]) for a/R > 0.08 at TV = 2400. The estimated 
Co increases proportionally with Cbd, as shown in Fig. 
ITifa): Cocr ^ 0.5Cbd- A small deviation (< O.Ol/a) from 
the line is almost independent of Cbd (see Fig. EJb)); 
therefore, it is considered a systematic error because the 
symmetric membrane at Cbd ~ has Co precisely equal 
to zero. 

The bending rigidity k is independent of Cbd- The re- 
sults coincide with the k estimated from the thermal un- 
dulations of the planar membranes (see Fig. [HTc)). The 
two methods have a slight dependence on Cbd with op- 
posite tendencies, but the dependences are smaller than 
the error bars. 

Alternatively, Cq can be estimated from the shape of a 
membrane strip. The membrane is connected by the pe- 
riodic boundary in one (x) direction, whereas it is open 
with edges in the other (y) direction (see the snapshot 
in Fig. [T3^a)). Because the membrane can freely bend 
in the y direction, the mean curvature should be Co. 
The flip of the orientation vector Ui of the membrane 
particles is not observed for the investigated parame- 
ters. Thus, the membranes can maintain the value of 
Co homogeneously even when the membrane has open 
edges. We calculated the membrane curvature using the 
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FIG. 14: (Color online) Parameter Chd dependence of (a), (b) 
the spontaneous curvature Co and (c) the bending rigidity k 
for the spin membrane model at e = 5. (a) The symbols (□) 
and (o) represent Co obtained from the membrane tube for 
fcbend = fctilt = 5 and fcbend = fctilt = 20, respectively. The 
solid line shows the line Cqct — 0.5Cbd. The snapshots of a 
membrane tube and strip are shown in the inset at A'^ = 1200 
and 400, respectively, for — 24a, fcbcnd ~ ^tiit = 20, and 
Cbd = 0.2. (b) The symbols (□, o) and (V, A) represent the 
relative spontaneous curvature Coo" — O.SCbd obtained from 
the membrane tubes and membrane strips, respectively, (c) 
The opened or closed circles represent the values estimated by 
the fits with two fit parameters (k, Co) or one parameter (k), 
respectively. The triangles (A, V) represent k estimated by 
the undulation analysis of the planar membranes with — 
1600 and 6400, respectively. 

second-order moving least-squares fit [5^1 with the weight 
function of the potentials Wcv{r) for strips of iV = 400 
with Lz/(J = 20 to 30. For a large bending rigidity 
K = 34:kBT (fctilt ~ 20), the resulting Co follows the 
relation Cq(j ~ 0.5Cbd better than that of the tube es- 
timation. However, it seems to underestimate Co for a 
small bending rigidity k ~ Qk^T (fctiit = 5) owing to 
large particle protrusions. We confirmed that the result- 
ing values were not sensitive to the shape of the weight 
in the mis fit. When a larger radius (rga = 2.5a and 
Tec = 5a) is used for the weight Wcv{r), the differences 
from the Co values calculated with the original weight 
are less than 5%. 

The bending elasticity generated by the bending and 
tilt potentials can be derived from the continuum theory 
as discussed in Ref. [2^. When the orientation vec- 



tors Ui are equal to the normal vectors of the membrane 
without tilt deformation, the bending and tilt energies 
are given by 

f/ev = JdA ^[(Ci - C^)2 + (C2 - C'o)^] 

+ ^{Cf + C',) (27) 

^ JdA <cnd + <lt (g^ ^^)2 

-Kc„d + <iit)CiC2 + C/o (28) 

in the continuum limit, where Ci and C2 are two prin- 
cipal curvatures of the membrane. The first and second 
terms in Eq. ((27)) are the contributions of the bending 
and tilt potentials, respectively. The spontaneous curva- 
ture of the bending potential is given by Cq = Cbd/j^nb- 
The nearest-neighbor distance fnb — 1.15a is obtained 
from the radial distribution function. By assuming a 
hexagonal packing of the molecules, the bending rigidi- 
ties generated by the bending and tilt potentials are esti- 
mated as Kbcnd/^BT = y3fcbcndWcv(nib) and K[-^JkBT = 
V5kti\tWcv{rnh)/2, respectively. The bending rigidity is 
given by their sum; i.e., k = k'^^^^^ + i^'aif Equation ([25)) 
gives the saddle-splay modulus k — —k and the spon- 
taneous curvature Cq = {K;,„„d/«end + '«iiit)}Cbd/r„b 

with Uo = (^Lnd + <ilt)(l/2 + <ilt/<e„d)Co- Thus, K 

and Co are estimated as 

K = (fcbcnd + 0.5fctiit)fcBr, (29) 

Co = {fc'bond/l-15o-(fcbcnd + fctilt/2)}Cbd, 

from Wcv(l-15o-) = 0.56. For fcbond = fctiit, Cq = 0.58Cbd- 
This relation explains the simulation results very well. 
The 16% overestimation of the factor (0.58) may be 
caused by the assumption of a regular hexagonal struc- 
ture for the fiuid state. 

Another method to estimate the spontaneous curva- 
ture Co in a simulation was proposed by Markvoort et 
al. [4^ . They made a sigmoidal shape of the membranes 
with two domains, which have opposite spontaneous cur- 
vatures. Then, Co was estimated from a comparison of 
the membrane shape with the energy minimum curve of 
the continuum theory (49l |. The accuracy of this method 
is likely to be similar to that of the curved strip since 
both the methods use the minimum energy shape of the 
membranes. 

Among these three methods, the Co estimation for 
tubular membranes can be applied even for a small 
bending rigidity k ^ lOk^T. It would be suitable for 
other solvent-free models, where the membrane can freely 
change the tube volume or the area difference between 
the two leafiets in a bilayer membrane. For membranes 
with explicit solvents or with slow flip-flop relaxation, the 
other two methods, with the curved membrane strip or 
sigmoidal membrane, would be easier to apply. 
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VI. SUMMARY 

We have investigated numerical methods for measur- 
ing the bending rigidity k and spontaneous curvature Cq 
of fluid membranes. For planar membranes, k, is esti- 
mated from the spectrum of the thermal undulations. It 
is found that estimate values show a large dependence 
on the upper-cutoff frequency gcut for the least-squares 
fits. Among the investigated fitting methods, the inverse 
power-spectrum fit with the extrapolation to gcut 
gives an accurate estimation. For tubular membranes, k 
is estimated from the stretching force and the spectrum 
of the thermal undulations. The estimated k, gives a rea- 
sonable agreement with the others for all of three meth- 
ods as well as for previous methods usiiig the anisotropic 
surface tension of a buckled membrane [39| and the ther- 
mal undulations of quasi-spherical vesicles [ssf . From a 
comparison of these methods, it is concluded that the 
inverse power-spectrum fit with the extrapolation is the 
best estimation method for simulations. 

The excess area Oox of tubular membranes is also in- 
vestigated. For short tubes, the calculated agrees 
with that obtained by Fournier et aUs perturbation the- 
ory [Tsj . However, with an increasing tube length, 
decreases in the simulation but increases in their calcu- 
lation. This difference may be caused by the finite-size 
effects in the simulations or the effects of the higher-order 
terms of the perturbation theory. 

The spontaneous curvature Co is measured from the 
axial force of tubular membranes and the average curva- 
ture of bent membrane strips. Both the methods provide 
a reasonable estimation. The methods investigated here 
are also suitable for other membrane simulation mod- 
els from atomistic or coarse-grained molecular models to 
large-scale meshless models. 
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Appendix A: membrane properties 

Here, we describe the parameter dependence of the 
properties of the tensionless membrane for the mis and 
spin membrane models. Figures [T5] and [T51 show the de- 
pendence of five quantities on curvature parameters (ka , 
fctiit, and fcbend) and attraction strength e, respectively. 
The membrane is in the fluid phase for all ranges of the 
parameters shown in the figures. 
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FIG. 15: (Color online) Curvature parameter dependence of 
(a) the intrinsic area Ao/Na'^, (b) area compression modulus 
Ka, (c) bending rigidity k, (d) line tension F, and (e) diffu- 
sion coefficient D for tensionless membranes at e = 4. The 
solid line with circles represents the data for the mis model. 
The dashed lines with squares and triangles represent data 
for the spin model at Cbd = for fcbcnd ~ ktnt and fcbond ~ 0, 
respectively. The solid and dashed lines in (c) show the linear 
fits for K/ksT = 2.3ka - 2.4 and K/ksT = 1.75A;tiit - 0.9, as 
well as for n/ksT — 0.52fctiit — 2.7, respectively. 



The intrinsic area Aq, the area compression modu- 
lus Ka, the line tension F of the membrane edge, and 
the particle diffusion coefficient D are almost indepen- 
dent of the curvature parameters when they are suffi- 
ciently large {ka > 10 and A;tiit > 15). The bending 
rigidity k is linearly dependent on the curvature param- 
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FIG. 16: (Color online) Parameter e dependence of (a) 
Ao/Na^, (b) Ka, (c) k, (d) F, and (e) D for the tensionless 
membranes. The solid line with circles represents the data 
for the mis model at ka — 10. The dashed line with squares 
represents the data for the spin model at ktut = fcbcnci = 20 
and Cbd ~ 0. The solid and dashed lines in (d) show the 
linear fits Ta/kaT = 1.15e - 0.2 and Ta/ksT = 1.12e - 0.6, 
respectively. 



cters (see Fig. \TE\i . Thus, k can be varied without 
changing the other membrane properties. For the spin 
model, the dependence of k on curvature parameters can 
be quasi-quantitatively explained by Eq. (|29p . derived 
from the continuum theory. The slope is only 17% or 4% 
higher than the theoretical prediction for fcbcnd = ^tiit 
and fcbcnd = 0, respectively. The line tension F linearly 
depends on e, whereas k is almost independent of e (see 
Fig. [T6|) . Thus, K and F can be separately varied by 
changing the potential parameters for the spin model as 
well as for the mis model. 



The intrinsic area Aq, the area compression modulus 
Ka — Aod-f/dA\A=Ao: and the diffusion coefficient D 
for a tensionless membrane are calculated from planar 
membranes using the method explained in Ref. [55| . The 
unit diffusion coefficient is Dq = jr^. The bending 
rigidity k is estimated using the extrapolation method 
for Eq. dUl) at TV = 1600. It is 10% higher than the 
values estimated in our previous paper |55| . where Eq. 
p6)) is fitted with (gcut/Tr)^ = 0.05. 



The line tension F of the membrane edge is calculated 
from the membrane strips with = 400, as follows: [67l - 



r = {{Pyy + PzzM'l - P^a,)LyLj2, (Al) 



since the length of the membrane edge is 2Lx and F is 
the energy per unit length of the membrane edge. The 
resulting F values coincide with the values estimated from 
the membrane pore in Ref. [55| . 



The difference between the values is less than 
ClfceT/cr. 
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